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Abstract Activator-inhibitor FitzHugh-Nagumo (FHN) equation is an example for reaction- 
diffusion equations with skew-gradient structure. We discretize the FHN equation using 
symmetric interior penalty discontinuous Galerkin (SIPG) method in space and average vec¬ 
tor field (AVF) method in time. The AVF method is a geometric integrator, i.e. it preserves 
the energy of the Hamiltonian systems and energy dissipation of the gradient systems. In 
this work, we show that the fully discrete energy of the FHN equation satisfies the mini¬ 
maximizer property of the continuous energy for the skew-gradient systems. We present 
numerical results with traveling fronts and pulses for one dimensional, two coupled FHN 
equations and three coupled FHN equations with one activator and two inhibitors in skew- 
gradient form. Turing patterns are computed for fully discretized two dimensional FHN 
equation in the form of spots and labyrinths. Because the computation of the Turing pat¬ 
terns is time consuming for different parameters, we applied model order reduction with 
the proper orthogonal decomposition (POD). The nonlinear term in the reduced equations is 
computed using the discrete empirical interpolation (DEIM) with SIPG discretization. Due 
to the local nature of the discontinuous Galerkin (DG) method, the nonlinear terms can be 
computed more efficiently than for the continuous finite elements. The reduced solutions 
are very close to the fully discretized ones. The efficiency and accuracy of the POD and 
POD-DEIM reduced solutions are shown for the labyrinth-like patterns. 
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1 Introduction 

The FitzHugh-Nagumo (FHN) equation was proposed for modeling the electrical impulses 
in a nerve axon as a simplified Hodgkin-Huxley model. Although the FHN equation is devel¬ 
oped as a model in physiology, it is used as a generic model that exhibits many phenomena 
in excitable or oscillatory chemical media. FHN equation is a singularly perturbed, param¬ 
eter dependent equation with a complex dynamics. In the literature the most known type of 
the FHN equation is the one which consists of a partial differential equation (PDE) for the 
activator and an ordinary differential equation (ODE) for the inhibitor. In this work we con¬ 
sider the so called diffusive FHN equation consisting of two or three PDEs with an activator, 
and one or two inhibitors. Depending on the relationship between the reaction and diffusion 
parameters, and the bistable nonlinearity, the specific patterns in one and two dimensional 
FHN equation are like traveling fronts, pulses, spots or labyrinth like patterns. Energy of 
the system plays an important role in the selection of these patterns. As a reaction diffusion 
equation FHN equation has been used to study the wave dynamics in excitable media and 
bifurcation phenomena. 

For large class of reaction-diffusion systems there exist Lyapunov functions, so that all 
the solutions converge to equilibria. But it is difficult to derive gradient structure for reaction- 
diffusion systems, because the reaction terms contain quite general nonlinearities. For cer¬ 
tain classes, gradient structures can be constructed. Among these are the activator-inhibitor 
type reaction-diffusion systems which are not order preserving and their linearized forms 
around the steady states are not self-adjoint. Yanigada ( [Yanagida] |2002b| introduced the 
concept of skew-gradient systems to investigate the stability of these systems. It was shown 
in ( [Yanagida] |2002b| under some restrictions on the parameters the diffusive FitzHugh 
Nagumo system and Gierer-Meinhardt system exhibit skew-gradient structure. In short, 
these activator-inhibitor equations consist of two gradient systems. 

In the following we consider the coupled reaction-diffusion system of the form 


Ti u t — d\Au + /(w,v) 
T 2 v t = d 2 Av + g(u,v) 


( 1 ) 


in a smooth bounded domain Q in W 1 (n < 2) and on a time period [0, T] with u = u(x,t) and 
v = v(x,t), x G £2, t £ [0, T]. Here T \, t 2 and d \, d 2 correspond to time scales and diffusion 
coefficients of u and v, respectively. The system Q has a skew-gradient structure if there 
holds ( |Yanagida] |2QQ2b| 



( 2 ) 


A particular example of the skew-gradient system is the diffusive FHN equation model¬ 
ing the transmission of electrical impulses in a nerve axon. For this equation, we can write 
/(m,v) and g(u,v) in equation 0 as a summation of separate functions of u and v as 

Ti« r = d\Au + f\{u) +f 2 (v) 
t 2 v, = d 2 Av + gi(u)+g 2 (v) 

where f\ ( u ) = u(u — /3)(1 — u), f 2 (v) = — v+ K, gi(u) = u, g 2 (v) = — yv + e and y > 0. The 
FHN equation 0 is an activator-inhibitor system, where u is the activator since it activates v, 
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i.e. leads to an increase of v. On the other hand, v is the inhibitor since it leads to a decrease 
in u and v. Moreover, FHN equation 0 is a skew-gradient system 0 through the energy 
functional 


E(u,v) = J f ^y|V(/| 2 - y |Vv | 2 + F(«,v) 


with the corresponding potential function 

u A (l+/3)w 3 
F(u,v) = --+ y 


j5u 2 yv 1 

'-—-uv+t— 
2 2 


ck 


-ev, 


(4) 


where V denotes the gradient operator. The first equation of 0 is a gradient flow with 
the potential F(w,v), and the second with the potential —F(w,v) ([Yanagida, j2002a). The 
energy functional E(u,v) does not correspond to the Lyapunov functional, because it is not 
necessarily non-increasing or non-decreasing in time, but the steady state («, v) = (0, iff) of 
a skew-gradient system is stable if it is a mini-maximizer of E(u,v) ([Yanagida, 2002b ). 

The stationary uniform solutions of ([ 3 ]) can be found by solving fi {<!>)+ fi{v) — 0 = 
gi (0) +g2(V / )- For one dimensional skew-gradient systems, the number of solutions of this 
equation determines the type of the solutions of the FHN equation 0. If it has only one 
intersection point, then it is called monostable case. However, if there is three intersec¬ 
tion points with two stable and one unstable fixed points, then it is called bistable case. 
For the monostable case, the solution of FHN equation ([ 3 } exhibits traveling pulses and 
for the bistable case, it exhibits traveling fronts. For two dimensional systems, under some 
conditions on the parameters, Turing patterns may exist, when the spatially homogeneous 
steady-state (w, v) is not a mini-maximizer of E(u, v) ^Yanagida, 2002b|. 

In this paper, we use symmetric interior penalty discontinuous Galerkin finite elements 
(SIPG) ( |Arnold| 1982[ Rivier e[|2008|) for space discretization and average vector field (AVF) 


method ( |Celledoni, E. and Grimm, V. and McLachlan, R. I. and McLaren, D. I. and O’Neale, 


|D. J. and Owren, B. and Quispel, G. R. W. 2012 Hairer and Lubich 2014| ) for time dis¬ 
cretization. Discontinuous Galerkin (DG) method uses discontinuous polynomials for the 
discretization of PDEs. DG methods can support high order local approximations that can 
vary non-uniformly over the mesh, by which fronts, pulses and layers can be captured bet¬ 
ter. The AVF method is second order in time and preserves for conservative systems, like 
Hamiltonian, the energy and for gradient systems the energy dissipation. We show that the 
AVF integrator combined with SIPG space discretization preserves the mini-maximizing 
property of the skew-gradient systems ( |Yanagida| |2002a| ) in the discrete form. 

FHN equation has front or pulse solutions in one dimensional domains depending on 
the parameters. We show in numerical simulations the existence of the multiple pulse and 
front solutions for the two and three component one dimensional FHN equations in the 
skew-gradient form. We also show for two the dimensional FHN equation in skew-gradient 
form the criteria of the existence of patterns like spots and labyrinths. The time evolution of 
the discrete energy shows that the energy remains constant in long term integration, which 
proves the stability of the solutions obtained by the numerical solutions. 

There has been significant development in the efficient implementation and analysis of 
the model order reduction (MOR) techniques for parametrized PDEs |Grepl| |2012| . The 
aim of the model order reduction is the construction of a reduced order model (ROM) of 
small dimension which essentially captures the dynamics of the full order model (FOM). 
The most known model order reduction technique is the proper orthogonal decomposition 
(POD). Even though the POD is a very successful MOR technique for linear problems, for 
nonlinear problems the computational complexity of the evaluation of the nonlinear reduced 
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model still depends on the dimension of the FOM. In order to reduce the computational cost 
several methods are developed so that the nonlinear function evaluations are independent of 
the dimension of the FOM and the computational complexity is proportional to the number 
of ROMs. Among the best known of them is the discrete empirical interpolation method 
(DEIM) ( [Chaturantabut and Sorensen j2010| which is adopted from the empirical interpola¬ 
tion method (EIM) (Barrault et al, 2004). The DEIM was originally developed for nonlinear 
functions which depend componentwise on single variables, arising from the finite differ¬ 
ence discretization of nonlinear PDEs. For the finite element discretization, the nonlinear 
functions depend on the mesh and on the polynomial degree of the finite elements. There¬ 
fore the efficiency of the POD-DEIM can be degraded. A modified version of POD-DEIM 
is developed in ( |Tiso and Rixen|[20l3] ) using the unassembled finite elements. This reduces 
the number of nonlinear function calls during the online computation, but the size of the 
nonlinear snapshots are enlarged, which increases the offline computational cost. In ( |Antil| 
|et al[ |2014| ), the assembled and unassembled finite element POD-DEIM are compared for 
parametrized steady state PDEs. In case of DG discretization, each component of the non¬ 
linear functions depends only on the local mesh in contrast to the continuous finite element 
discretization where the nonlinear function depends on multiple components of the state 
vector. Therefore the number of POD-DEIM function evaluations for DG discretization is 
comparable with the finite difference discretization. 

The paper is organized as follows. In Section[2] the discretization of FHN equation ^ in 
space by SIPG is given. In Section]?] the fully discrete formulation of the FHN equation ([3| 
is presented using the gradient stable AVF time integrator. The Energy analysis of the fully 
discrete scheme is given in Section]?] In Section]?] numerical simulations are shown for the 
FHN equation exhibiting front and pulse solutions in one dimensional domains. We present 
in Section [6] the numerical results for the three component FHN equation in skew-gradient 
form with multiple pulse solutions in one dimensional domains. Turing patterns for two- 
dimensional FHN equation are shown in Section [ 7 ] We show in Section [ 8 ] the effectiveness 
and accuracy of the reduced order models for SIPG discretization with application of the 
POD-DEIM to the nonlinear term. The paper ends with some conclusions. 


2 Space discretization by discontinuous Galerkin method 

In last twenty years, the DG methods gained an increasing importance for an efficient and 
accurate solution of PDEs. Although continuous finite elements (FEM) require continuity of 
the solution along element interfaces, DG does not require continuity of the solution along 
the interfaces. In contrast to the stabilized continuous Galerkin finite element methods, DG 
methods produce stable discretization without the need for extra stabilization strategies, 
and damp the unphysical oscillations. Due to the local structure of DG discretization, DG 
methods are parallelizable and adaptive meshing techniques can be implemented efficiently. 
The DG methods combines the best properties of the finite volume and continuous finite 
elements methods. 

In this section, we will describe the DG discretized semi-discrete (continuous in time) 
form of the FHN equation © with homogeneous Neumann (zero-flux) boundary condition. 
The classical (continuous) weak solution of ^ solves for t £ (0, T] the variational formula¬ 
tion 


(TiM;,wi)+a(di;M,wi)-(/i(M),wi)-(/ 2 (v),wi) = 0, Vwi, 

(t '2Vt,w 2 ) + a(d 2 -,v,W2)-(gi{u),W2)-(g2(v),W2 ) =0, Vw 2 


( 5 ) 
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with the initial conditions satisfying 


0(0),wi) = Oo,wi), (v(0),w 2 ) = (vo,w 2 ), 


where (•,•) := (-j-)n is the L 2 inner product over the domain Q, w\ and w 2 are the test 
functions and a(d\u,w) — (dVu,Vw)n stands for the classical bilinear form. 

Let Sh be the disjoint partition of the domain Q with elements (triangles) G £h, 

where N e [ is the number of elements in the partition. On £/*, we set the discrete solution and 
test function space 


Dk — D k (£ h ) { w ^ L 2 (Q) : we G ^k(E) \/E G £h }, 

where P^(£) is the space of polynomials of degree at most k on E G Multiplying by 
w\ and W 2 and integrating by using Green’s theorem over each mesh element, we obtain the 
following semi-discrete variational formulation: \/t G (0, T], find Uh(t ), Vh(t) G Dk satisfying 

( T i^r’ lv i) + a h{d\\u h ,W]) - (fi(u h ),wi) - (f 2 (v h ),wi) = 0, Vwi G Dk, 

} dv \ ( 6 ) 

+ a h (d 2 \v h ,w 2 )-(g\(u h ),w 2 )-(g 2 (y h ),w 2 ) 0, Vwi G D k , 


where the SIPG bilinear form ah(d ; w, w) is given by 


a h (d-,u,w)= Y / 
Eee h 


dWu • Vw — 


£ f{dVu}-[w\ds 

6 

Y f {dVw}-[u]ds+ Y /{Vw}-[w]ds, 

n e Je 


eer" 


where h e denotes the length of an edge e , L^° denotes the set of inter-element faces (edges), 
[ ] and { } stand for the jump and average operators respectively. Introducing the degrees of 
freedom N \—Ni oc x N e i, where Ni oc denotes the local dimension on each element depending 
on the polynomial order k , the semi-discrete DG solutions of © are of the form 

N N 

Uh{t ) £wi(fM > v /*(0 = ( 7 ) 

z=l i=l 

where u (t) := (u\(t),. . .,u^(t)) T and v(t) := (vi (t),. .. ,v^(t)) r are the vectors of time de¬ 
pendent unknown coefficients of u e and v/*, respectively, and (j) := (</>i,..., 0^v) r is the vector 
of basis functions. Plugging ([ 7 } into the scheme ([6} and choosing = w 2 = <t>i, i — 1, • • • ,N, 
we get the system of 2 x N dimensional ODEs for the unknown vectors u and v as 


TiMu t + S u u - F\ (u) - F 2 (v) = 0, 
T 2 M\ t + S v \ - G\ (u) - G 2 (u) = 0, 


where M, S u , S v G R NxN are the mass matrix and stiffness matrices, respectively, and the 
remaining are the vectors in W N of the unknowns u and v. 
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3 Time discretization by the average vector field method 

Energy stable time discretization methods preserve the dissipative structure of the numerical 
solution of gradient flow equations and skew-gradient systems like the FHN equation. The 
small values of the diffusion parameters leads to stiff systems after spatial discretizations. 
Implicit/implicit-explicit methods are developed since the explicit methods are not suitable 
for stiff systems and the fully implicit systems require solution of non-linear equations at 
each time step. In the semi-implicit schemes, the linear stiff part is treated implicitly and 
the nonlinear part explicitly, so that at each time step a linear system of equations is solved. 
However, these methods do not preserve the energy dissipation of the system. Implicit Eu¬ 
ler method and average vector field (AVF) method are energy stable time discretization 
techniques which are robust with small diffusion parameters. The AVF method is the only 
second order implicit energy stable method and it preserves energy decreasing property for 
the gradient systems and for systems with Fyapunov functionals like the FHN equation. In 
this work, we apply AVF method to solve the system of ordinary differential equations 0 
arising from the semi-discretization of the model problem 0- 

We split the time interval [0, T] into J equally-length subintervals {tk ~\, h\ with 0 = to < 
t\ < ... <tj — T with the uniform step-size At = h~ 4-i, k — 1,2,...,/. The AVF method 
for an arbitrary ODE y = f(y) is given by 

+ (9) 

where y n ~y(t n ) for n— 1,...,/ and yo = y(/o)- Applying the AVF formulation to the 
system of ODEs 0, in space SIPG discretized fully discrete formulation of the equation 0 
reads as 

y5„j u n + At /F l( ^i + (l-«K)d« 

+ At J F 2 (£v„ + i + (l-£)v„)d£, 

° , ( 10 ) 
(yiM - y S v ^j v„ + At JGi (|u„+i + (1 -1 )u„)d£ 

1 

+ At J G 2 (^v n+ i + (1 —%)\ n )dt;. 
o 

The fully discrete system of nonlinear equations © is solved by Newton’s method on each 
time-interval (t n -\,t n \ , n = 0, 1,..., / — 1 . 

4 Energy analysis of fully discrete scheme 

It was proved in ( |Yanagida| |2QQ2a| that a steady state of a skew-gradient system is stable if 
and only if it is a mini-maximizer of the energy functional {12} : 

- If u = ^ is a local minimizer of E(u,rj) and v = rj is a local maximizer of E(E ,, v), then 
(w,v) = (£,rj) is a mini-maximizer of E{u,v). 


At . 

Z\M + ~^S U ) u n+ i 


At 

t 2 M+—S v j v„+i 
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- If E(!,v) <£(!,!?) < E(u, tj) for any neighborhoods w, v of § and rj, resp., then we 
say that (w, v) = (£, 77 ) is a mini-maximizer of is (w, v). 

We will examine the relation between the stability of a steady state solution (£, 77 ) of 0 
and the mini-maximizing property of the critical point of the energy functional for the 
SIPG-AVF fully discrete system. 

Let u n and v n denote the solutions of 0 at t — t n . Steady state solutions (u n ,v n ) m 
(4(x),v(x)) Of0 satisfies 


d\AE, +/i(^) + fi{v}) — 0, 
J 2 zir}+gi(^)+g2(?7) = 0. 


( 11 ) 


Then the solution of GD is a critical point of the following energy functional 

E(u n ,V n ) tss — II Vm w || L 2 ( £/j ) — ll^ v n|lL 2 (£/j) + {E{ u m V n ), 1) 

(J(/i 


( 12 ) 


T V ( [w w ]) e + ([ M n]j[ w n]) 

v 2/!e 

where F(-, •) is the potential function 0 . Fully discrete variational formulation reads 
1 l 

^ t ( U n +1 -Mn>Wt) = - + M n ,Wl)+ J (f\(£,U n+ \ + (1 - § )u n ), Wi )(/§ 

0 

1 

+ f (/2(^V n+ i + (1 — 

0 , (13) 

^2 1 f 

^(V n +I -Vn,W 2 ) ~ 2 fl /i(^ 2 ;Vn+l +V„,W 2 )+ / (gl (§ M„+i + (1 - £ )u n ), W 2 )^ 

0 

1 

+ J(g2(%Vn + l + (l-Z)Vn),W 2 )dZ 
0 

By taking W] w n+ i - and w 2 v n+ i - in and using the identity (a + b,a — b) = 
(a 1 — b 2 , 1) and the bilinearity of < 2 ^, we get 

1 5 1 Mn) ~ ^ (d\ \ U n -\- 1 , U n + 1) T ~CLfo {d\ \ U n , U n ) 


+ (/l (%U n+ i + (1 -E,)u n ),u n+ 1 - 

+ (/2(^V n+ l +(1 -^)v n ),M n+ l -Un)d£,, 

T 2 11 

^(v«+i -v„,v n+ i — v„) = --a h (d 2 ',Vn+i,v n+ i) + -a h (d 2 -,v n ,v n ) 

+ {gl(Z u n+l + (1 -|)«n),V„ + l -v n )d£, 
+ {g2{&n+i + (1 -|)v„),v n+ i -v n )dE,. 


(14) 
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Using the Taylor expansion for the terms /i, / 2 , g 2 in and subtracting F(u n ,v n ) 
from F(w n+ i, v n+ i), we get 


dF 

F(u n +uv n +i)-F(u n ,v n ) « —(§u n +\ + (1 -£)M„,v n )(w n+ i -m„) 

<7W 

dF 

H ^ (<■? ^«+l T (1 <5 v n ) (v n _)_i v w ). 

<7V 


(15) 


Using the derivatives ^ = /i(w) +/ 2 (v), ^ = — (gi(u) + g2(v))> which are the skew- 
gradient conditions and plugging the identity G3 into the system we obtain the 
relation 


E(u n + l,V n +i) E(u n ,V n ) ~ ^ || M n+l u n\\ L 2 (n)+ ^ || v «+l v «IIl 2 (I2) 




(16) 


Now, at the steady state point (u„,v„) = ||v„+i - = 0 if v is 


fixed to 


7] (x), leading according to the relation ( fib} to £(<^+i ,7]) < is (£ n , 7]). In a similar way, if u 
is fixed to £ (x), we obtain is (£, 7] n+ i) > is(£, 7] n ). Thus, the first equation of 0 describes 
a gradient flow with E(u,rj) and the second equation of ^ describes a gradient flow with 
—is(£,v), meaning that (u n ,v n ) = (£,tj) is stable as a steady state of 0 since it is a mini¬ 
maximizer of is(w, v) (Yanagida, 2002a| . 


5 Traveling fronts and pulses of FitzHugh-Nagumo equation 


Localized structures like fronts and pulses are most-well known one-dimensional waves in 
reaction diffusion systems ( [Chen and Hu[ |2014| . Fronts connect two different state of a 
reaction-diffusion system with a bistable nonlinearity. Pulses exist far away from the homo¬ 
geneous equilibrium and results from the balance between the dissipation and nonlinearity. 
Existence of fronts and pulses for the FHN equation 0 are shown in ( [Chen and Hu||2014} . 
In this section we derive conditions for the parameters of the FHN equation ^ exhibiting 
traveling fronts or traveling pulses which depend on mono/bi-stability of the homogeneous 
steady state solutions of 0 

f(u) — v = 0, yu — <5v + e = 0. (17) 

Here, we consider the bistable cubic nonlinearity f(u) — u(u — j3)(l — w), with 0 < /3 < 
Eliminating v = from {IT} and fixing the parameters /3 = ^, 7=1, £=^,we obtain 


ic’ — —— u -\- { —— -j— ) u -\- 


25 


25 


lOy 


= 0 . 


(18) 


For the bistable case, the equation ( [18] ) must have three distinct reel roots, or equivalently 
its derivative equation 


3 u 2 


54 (2 

25“ + \25 


(19) 


must have two distinct roots. In other words, we need the condition 


A = 



2 

-4-3- 


2 1 
25 + y 


2316 

“625“ 



Hence, for y > 3.2383, the equation 0 attains the bistable case, as a result, the solution of 
0 is traveling fronts. Otherwise it will be monostable and the traveling pulses will occur. 
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5.1 Traveling fronts 

We consider the FHN equation 0 on Q = [—60,60] with the parameters d\ = 1, d 2 = 1.25, 
T i=s 12.5, /5 = y = 8 and £ = 0.7. We set the initial conditions u(x, 0) = tanh(x), v(x, 0) « 
1 — tanh(v), and the spatial mesh size Ax = 0.1 and the temporal step size At — 0.5. The 
traveling front solutions can be clearly seen from Fig. [T] whereas the discrete energy is 
decaying very slowly. 


t=100 t=250 



x x t 


Fig. 1: Traveling front solutions and evolution of the discrete energy for Example 


5.1 


5.2 Traveling pulses 


All parameters other than y= 0.8 are the same as for the traveling front solutions in Example 
5.1 Initial conditions are taken as u(x, 0) = tanh(x) and v(x,0) = —0.6. The traveling pulse 
solutions obtained for Ax = 0.1 and At — 0.1 are shown in Fig. [2] The discrete energy 
remains constant after an oscillation. 


t=20 t=60 



x x t 


Fig. 2: Traveling pulse solutions and evolution of the discrete energy for Example 


5.2 
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6 Three component FitzHugh-Nagumo system 


The rich dynamics of localized structures like multi-pulse and multi-front solutions were 
recently investigated for singularly perturbed three component reaction-diffusion systems 
( |Or-Guil et al||1998||van Heijster et al||2008||van Heijster and Sandstede||201 1) 

u t u xx + u-u 3 -£(av + /3s + y), 

1 

tv,= ^v« + k-v, (20) 

d 2 

0s t = -^s xx + u-s, 

£ l 

where 0 < £ 1, t, 0 > 0, d > l, and a and j 3 denoting the reacting rates and the constant 

y is the source term. The system <(20|» was originally introduced as a model for gas discharge 
dynamics ( [Or-Guil et al[|1998) . It became a standard model to study the dynamics and inter¬ 
actions of spatially localized structures like pulses and fronts in one dimensional reaction- 
diffusion equations. The three component model ( |20) is also considered as augmented FHN 
equation with a second inhibitor w which diffuses more rapidly than the second inhibitor v. 

The last two equations of ([20) imply that all the components of the stationary solutions 
are equal. Elimination of the components v and s leads to the equation 

u 3 — u £(ocu f3u = 0 , ( 21 ) 


which has three different roots if its derivative equation 3u 2 — l + e(a + /3) has two distinct 
roots, i.e. £(a + /3) < 1. In this case, it is called bistable medium since two of the roots of 
( [21) is stable. Hence, ( [20) has front solutions, and depending on the initial conditions we 
can get one-front, two-front or multi-front solutions. 

On the other hand, the three component FHN equation ( |20) has skew-gradient structure 
if there hold 


( 22 ) 


£oc _ 1 e/3 _ 1 
t 0’ t O' 

In addition to the skew-gradient condition ( [22) , if e(a + /3) > 1 is satisfied, then we obtain 
one stationary point. In this situation, the stationary point is always unstable since the Jaco¬ 
bian matrix around this stationary point of ( [20) satisfies det(/) < 0 and oscillatory traveling 
waves exists satisfying the skew-gradient condition ( [22) . 

In the following numerical examples we present one pulse (two-front), two pulse (four 
front) and multi-pulse (multi-front) solutions of ( [20) similar to those in ( [van Heijster et al[ 
[2008) for different set of parameters satisfying the skew-gradient condition ( [22) under the 
homogeneous Neumann boundary conditions on both ends. In all examples we have taken 
£ = 0 . 01 . 


6.1 One pulse solutions 


We consider the space-time domain x G [—1000,1000] and t G [0,200] with Ax — At — 0.5, 
and the initial conditions are given by 

x G [-50,50] 

— 1 otherwise 


u{x, 0) = | (j 


v(x,0) — s(x,0) — 0. 


We use the set of parameters (a,/3, y,d, t, 0) = (3,1,-0.25,5,100/3,100). The solution 
profiles at t = 25 and at the final time t — 200, together with the energy plot, are shown in 
Fig H] 
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t=25 t=200 Energy 



x x t 


Fig. 3: One pulse solutions and evolution of discrete energy for Example 


6.1 


6.2 Two pulse solutions 


For the two pulse solutions, we use the same settings as in Example [67T| with the final time 
T — 100. We give the related results in Fig.[4]for the initial conditions 


f 1 , if* G [-350,-150] U [150,350], 
\ — 1 , otherwise. 


v(*,0) — s(x,0) — 0. 


t=20 t=100 Energy 




Fig. 4: Two pulse solutions and evolution of discrete energy for Example 


6.2 


6.3 Multi pulse solutions 

Interaction of multi-fronts are shown in Fig. [5] for the set of parameters (a, /3, 7 , d, T, 0) m 
(100,100, —0.25,5,1,1), and for the initial conditions 


«M) = { _? I otherwise. 501 * ’ v{xJ)) = s{xJ)) = °' 
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t=10 t=50 Energy 



x x t 


Fig. 5: Multi pulse solutions and evolution of discrete energy for Example 


6.3 


7 Pattern formation for the FitzHugh-Nagumo equation 


In this section, we will study Turing’s pattern formation of FHN equation ([3]). In order to 
determine what kind of spatial patterns is formed for Turing systems, an analytic criterion 
is developed in (Marquez-Lago and Padilla, 2014). Here, we will give this result and Turing 
instability conditions for the FHN equation ([3) with two different sets of f\ (u), / 2 (v), gi(u), 
g 2 (v) in two-dimensional spatial domain Q. Alan Turing proposed this relationship between 
the parameters of reaction-diffusion equations and pattern formation ( [Turing) |1952| . He 
showed that under certain conditions any stable steady state can be converted into an un¬ 
stable state. This is known as diffusion driven instability or Turing instability. The Turing 
instability conditions for the system are given by 


d_h 

du 


dfi dg2 
du dv 
dg2 _ df2 dg]_ 
dv dv du 

<h%- +d <¥ 

du dv 

d ^+“'¥i 

du dv J 


< 0 , 

> 0 , 

> 0 , 

> Ad\d2 



dg2 

dv 


dfl dgi \ 

dv ‘ du J 


(23) 


which have to be evaluated at the steady state solutions. The first two conditions of ( |23) are 
already the stability conditions of the steady state solutions of The last two conditions 
make this stable steady state solutions unstable under the presence of the diffusion terms. 

Turing patterns for FHN equation occur in form of spots or labyrinth-like patterns, de¬ 
pending on the number of global minimums of the Lyapunov energy functional (see 
(Marquez-Lago and Padilla] [2014) )). When the Lyapunov energy functional (|4]) of the sys¬ 
tem has a unique global minimum, then patterns appear as spots, otherwise, labyrinth-like 
patterns occur.. 


7.1 Pattern examples 

We fix the parameters T\ = T 2 = \,d\= 0.00028, and we take the bi-stable nonlinear term 
/1 ( u ) = u — w 3 , fiiv) — K — v,g\(u) —u and gi(y) — — v. First of all, note that the reaction 
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terms include an even term K. So for K — 0, the pattern formation will be labyrinth-like 
patterns, and for nonzero K the results will be spots. We have to determine the values of d 2 
so that the Turing instability conditions ( [23] ) are satisfied. The only steady state solution of 
with the given conditions in this problem is (uo, vo) ~ (—0.368403, —0.368403). When 
we evaluate the first order partial derivatives at this point, we get 


d_h 

du 


0.592838, 


dfi 

dv 


- 1 , 



dg2 

dv 


The first two conditions of ( [23] ) are clearly satisfied. In order to satisfy the third and last 
conditions of ( [23] ), we need d 2 > 0.000472 and d 2 > 0.002242, respectively. Hence, for 
Turing instability, we obtain d 2 > 0.002242 . 

We take d 2 — 0.005 and consider the domain (jt,y) E [—1,1] x [—1,1] and t E [0,200] 
with the uniformly distributed random number between — 1 and 1. Space mesh size and time 
step size are taken as Ax — 1/16 and At — 1/10, respectively. The steady state is reached 
at t — 200. The corresponding solution profiles and energy plots of spot and labyrinth-like 
patterns are given in Fig. [6] and Fig. [7] respectively. 



Fig. 6: Spots: Solution profile of the first component u (left), second component v (middle) 
and evolution of discrete energy (right) for Example [771] 



Fig. 7: Labyrinth-like patterns: Solution profile of the first component u (left), second com¬ 
ponent v (middle) and evolution of discrete energy (right) for Example [771] 
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8 Reduced order model approximation for the FitzHugh-Nagumo system 


For MOR we consider the following special form of the FHN equation 

u t — d\Au + f{u) —v+K 
v t — d2Av — v-\-u 

The SIPG semi-discretization of([24]) leads to the following system of ODEs 

Mu t — S u u + F(u) — Mx + K 
Mx t S v x - Mx + Mu 


(24) 


(25) 


with the stiffness matrices S U ,S V E R NxN , M E R^ Vx ^ v is the mass matrix, K E R N is the 
constant vector related to the integral containing the parameter K and F E R N is the vector 
of the bi-stable nonlinear terms of u = (u \, U 2 ,..., un). The number N = Ni oc x N e [ stands 
for the degrees of freedom in DG methods, where N[ oc and N e i represent the number of 
local dimension in each element (triangle) and the number of elements, respectively. The 
solutions of ( [25] ) are of the form 

N N 

u(t) = Y Ui(t)<j>i(x) = <j>u(t) , v(t) = Y W = <M0 (26) 

i=l i= 1 

where E R are the unknown coefficients and fa are the DG basis functions. 

Because the computation of the Turing patterns for ( [25] ) is very time consuming, we 
construct the reduced order model (ROM) of small dimension ( [25] ) by utilizing the POD 
method jKunisch and Yolkwein[|2Q0l] ). 


8.1 Reduced order model 

The ROM for ( [25] ) of dimension k <^N is formed by approximating the solutions u(t) and 
v(t) of the FOM ( [25] ) in a subspaces spanned by a set of M-orthogonal basis functions 
{WuALi an d { YvJu=i of dimension k in R N and then projecting onto that subspace. The 
approximate ROM solutions have the form 

k k 

u{t) « £w;0)w > v (0 ~ (27) 

i— 1 i=l 

where u (t) = .. ,%(t)) r and x(t) — (vi(f),... ,v&(0) r are the coefficient vectors of 

the ROM solutions. The reduced basis functions and {l J/ V: i} of the ROM solutions in 

dz) are constructed with DG finite element basis functions {fa}^ 

N N 

= W = Y 'KjAi*) = (28) 

7=1 7=1 

where the coefficient vectors of the i — th reduced basis functions \f/ u j and y/ V j are located at 
the i — th columns of the matrices % = [^,.,1,..., E R Nxk and % = ,..., E 

R Nxk , respectively. 

The M-orthogonal reduced basis functions {tyuj} and {i// VjZ }, i — 1,2,... ,&, are com¬ 
puted by the POD ( [Kunisch and Yolkwein[|2Q0l| . We consider the snapshot matrices °7/ — 
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[u 1 ,..., u J ] and Y = [v 1 ,..., \ J ] in R NxJ , where each member of the snapshot matrices 
and *¥ are the corresponding coefficient vectors of the discrete solutions {u l } J i=l and {v l } J i=v 
respectively, of the FOM ( [25] ) at the time U,i = 0,1,..., /, with u l w ufc) and v l w v(fj). Then, 
for w G {w, v}, the M-orthogonal reduced basis functions { 1 / 4 ^}, i = 1,2,..., k, are given by 
the solution of the minimization problem 


min 

Vw,l,---iVw,k 



k 


Vw,i)l2(i2)¥w,i 


2 

L 2 (C2) 


subject to (VAi',;: Vwj)l 2 (q) — f j — k, 


where 5,- ; is the Kronecker delta. The above minimization problem is equivalent to the eigen¬ 
value problem 


W t M¥ u ,j = 


yy T M%,,i = , i = l,2,...,k 


(29) 


for the coefficient vectors and of the POD basis functions and \ff v j, respec¬ 
tively. Defining — Rtf/ and V — RV (R T is the Cholesky factor of the mass matrix M), 
we obtain the equivalent formulation of ( [29] ) as 


, vy 1 %„i = , i=l,2,...,k 


4?T 


(30) 


where v — R x R., ,i. The solutions of ( [30} are obtained as the first k left singular vec¬ 
tors 'Fu .j = ^ U i and %.i = in the singular value decomposition (SVD) of 9/ and 
respectively, as 

& = SuZufi, y = ^Pj, 

where the diagonal matrices E u and E v contain the singular values o U j and a v j on the di¬ 
agonals, respectively. Using the coefficient vectors of the POD basis 

functions are computed as 


11/ _d-1iia 




i— 1 , 2 ,.. .k. 


In addition, using the expansions ([26}, © and (|28]>, we have 


u = %U , V = 


(31) 


For the construction of the ^-dimensional ROM, we substitute the relations into 
the system ( [25} and we project onto the reduced spaces spanned by { 1 / 4 , 1 ,..., 1 / 4 ^} and 
{ 1 / 4 , 1 ,..., 1 / 4 ^}, respectively, leading to the system 


U, = Au + Tf F(%u) - M u \ + ff K 

\ t = £v v + M v u 


with the reduced matrices 

A = y,jA%, b = % t b%, m u = vJm%, m v = % t m%. 
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8.2 Discrete empirical interpolation method (DEIM) 


Although the dimension of the reduced system ( [32] ) is small, k<^N, the computation of the 
nonlinear term 


N(u) = 'FJf(%u) 


(33) 


still depends on the dimension N of the full system. We apply the DEIM (Chaturantabut and 
Sorensen 2010) to approximate the nonlinear function F(^u) from a subspace generated 
by the non-linear functions. Let & — [iq ,7^,... ,Fj\ £ R NxJ denotes the snapshot matrix of 
the nonlinear functions at each time instances computed as the solution of the full 

system ( [25] ) with Ft = F( x F u u(ti)), i— 1,2,... ,7. Using the SVD of the matrix we can 
find m<^N orthogonal basis functions {Wj}^ spanning the m -dimensional subspace of &. 
Then, with W = [W \,..., W m ] £ R Nxm , we can approximate the nonlinear function by 




(34) 


with the corresponding coefficient vector s(t). We note that the system © is overdeter¬ 
mined. Thus, to compute the coefficient vector s(t), we take m distinguished rows from the 
system Ws(t) through the projection using a permutation matrix P — \e ^ x ,..., ep m ] £ R Nxm 
where e^, is the i -th column o f the identity matrix I £ R A/xA ^ and computed by Algorithm [l] 
(Chaturantabut and Sorensen 2010) so that P T W is non-singular. 


Algorithm 1 DEIM Algorithm 

Input: POD basis functions {W/}^ =1 



Output: Index Vector p = \p\ ,..., p m \ 

T and Permutation Matrix P = [e ^ x ,..., ep m \ 

1 

INPUT: {w;}” , C R N 



2 

OUTPUT: p= [pi,...,p m ] T eR m , P eR Nxm 

3 

[\p\,Pi] = max{|Wi|} 



4 

W = [W l \,P=[e Pl ],p=[/p l \ 



5 

for i = 2 to m do 



6 

Solve (P r W)c = P T Wj for c 



7 

r = Wi - Wc 



8 

[\PlPi] = max{|r|} 



9 

W^[WWi],P^[Pe Pi ],p+- 

>" 

Pi. 


10: end for 




The projection of ( [34] ) leads to the system 

s{t) = (P T W)~ l P T F(%u(t)). (35) 

Then, using ( [34] ) and ( [35] ), the non-linear term ( [33] ) can be approximated as 

N(u(t))*N(u(t)) = QF, (36) 

where the matrix Q = WjW ( P T W)~ l £ R^ xm is precomputable and F = P T F('¥ u u(t )) £ IR m 
is the m-dimensional non-linear vector which can be computed in an efficient way. 

The DG requires only computation of the integrals on a single triangular element, which 
is not the case in continuous finite elements where all the interior degrees of freedoms are 
shared by usually 6 triangular elements (see Fig. [8]). The unassembled finite element ap¬ 
proach is used in ( [Antil et al||2014| , so that each DEIM point is related to one element. This 
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reduces the online computational cost, but increases the number of snapshots and therefore 
the cost of the offline computation. Due to its’ local nature, the DG discretization is automat¬ 
ically in the unassembled form and it does not require computation of additional snapshots 
(see Fig. [8). 




Fig. 8: Connectivity of degrees of freedoms for linear basis 


For the solution of the nonlinear system by Newton’s method, the entries of the Jacobian 
Jf £ M. NxN of the non-linear terms are of the form for /, j — 



Because the DG basis functions are defined only a single element and they vanish out¬ 
side that element, the integral terms of the Jacobian matrix (Jf)ij vanish on the triangular 
elements E t{ where the basis function (j)j is not defined. Unlike the continuous finite elements 
where the Jacobian matrix contains overlapping blocks, Jacobian matrix in DG appears in 
block diagonal form, and has the form 

= vJJf'Pu , ^V(u) = Q(P T J F )%. 

We note that (P T Jf) £ R mxA4 i s the matrix whose z-th row is the pi- th row of the Jacobian 
Jf, i = 1 , 2, ..., m, and in each row of the Jacobian there are only Ni oc nonzero terms because 
of the local structure of the DG. 


8.3 Numerical results 

We consider the FHN system ([24} with labyrinth-like patterns on Q = [—1, l] 2 with zero 
flux boundary conditions and with the initial conditions as uniformly distributed random 
numbers between —1 and 1. We set the parameters d\ — 0.00028, d2 = 0.005 and K = 0. 

In Fig. [9] the decay of the singular values related to U , V and the nonlinear term F are 
shown. In Figs. [To] [Tl] [12] we demonstrate that the ROM solutions and energy plots for 
k — 1 POD and m — 10 DEIM basis functions are almost the same as the FOM solutions in 
Fig[7] The computational efficiency of the POD-DEIM is obvious from Fig.[l3| 
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Fig. 9: Decay of the singular values related to U, V and nonlinearity F 



7 POD, 10DEIM 



-1 -0.5 0 0.5 1 


Fig. 10: Solution profiles of the component u at the steady state 



Fig. 11: Solution profiles of the component v at the steady state 


We demonstrate the computational efficiency of the POD-DEIM algorithm for linear 
and quadratic DG elements. In Table[l]we give the CPU time and the speed-up factors Spod 
and Sdeim of POD and POD-DEIM, respectively, 

CPU time for FOM CPU time for FOM 

P0D ~ CPU time for POD ROM ’ DEIM ~ CPU time for POD-DEIM ROM' 

In all computation only one Newton iteration is needed. The results show clearly that 
as the mesh size increases, the speed-up factors increases, which shows the efficiency of 
POD-DEIM method. In Table [ 2 ] we show the computed mean relative L 2 errors between 
FOM and POD, and between FOM and POD-DEIM solutions which are acceptable because 
POD-DEIM is an approximation of the non-linearity in contrast to the POD. 
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Full 7 POD 7 POD, 10 DEIM 



0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 

t t t 


Fig. 12: Energy plots 



Fig. 13: CPU times 


Table 1: The computation time (in sec) and speed-up factors for k = 10 POD basis functions 
and m — 50 DEIM basis functions on different grids and with different polynomial degrees 


Polynomial degree 


P= 1 



P = 2 


Mesh number 

1 

2 

3 

1 

2 

3 

# triangles 

128 

512 

2048 

128 

512 

2048 

# nodes 

81 

289 

1089 

81 

289 

1089 

# DoFs 

384 

1536 

6144 

768 

3072 

12288 

Full 

9.92 

26.25 

124.37 

29.42 

81.53 

423.97 

POD 

7.65 

13.50 

33.64 

22.23 

38.67 

116.82 

POD-DEIM 

7.11 

9.64 

12.67 

16.29 

19.12 

31.23 

SpOD 

1.34 

2.00 

3.83 

1.34 

2.15 

3.71 

Sdeim 

1.40 

2.90 

9.81 

1.86 

5.50 

15.39 


9 Conclusions 

Neurons, the building blocks of the central nervous system, are highly complex dynami¬ 
cal systems. The FHN equation as the simplified version of the Hodgkin-Huxley equation 
models in a detailed manner activation and deactivation dynamics of a spiking neurons. The 
FHN equation can describe the bifurcations with the variation of the key parameters for neu¬ 
ron dynamics. Within time FHN equations became a favorite model for simulation of wave 
propagation in excitable media, such as heart tissue or nerve fiber. Understanding the com¬ 
plex behavior of patterns in neuroscience can have big impact for dealing and preventing 
with various diseases. In this paper we have shown that patterns like pulses, fronts, spots 
and labyrinths for the FHN equation are computed accurately using the structure preserving 
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Table 2: Mean relative L 2 errors between the solutions of the full system and the ones for 
POD system and POD-DEIM system with k — 10 POD basis functions and m — 50 DEIM 
basis functions 


Polynomial degree 


P= 1 



P = 2 


Mesh number 

1 

2 

3 

1 

2 

3 

POD (Comp, u) 

3.8e-2 

1.2e-2 

1.0e-2 

8.9e-3 

l.le-2 

l.le-2 

POD-DEIM (Comp, u) 

3.6e-2 

1.5e-2 

1.5e-2 

1.2e-2 

2.5e-2 

2.9e-2 

POD (Comp, v) 

1.4e-2 

5.6e-3 

4.7e-3 

6.0e-3 

6.1e-3 

6.4e-3 

POD-DEIM (Comp, v) 

1.4e-2 

6.0e-3 

5.4e-3 

6.6e-3 

8.1e-3 

9.1e-3 


time integrator AVF in combination with DG finite elements in space. MOR can dramati¬ 
cally reduce simulation costs by preserving the behavior of parametrized PDEs, which was 
demonstrated here for the FHN equation with Turing patterns. In the literature MOR is ap¬ 
plied usually in connection with the finite differences and continuous finite elements. We 
have shown that the DG discretization in space can save the computational cost and due its 
local structure DG is more efficient than the continuous finite elements. 

Quantitative methods are indispensable to structure the clinical, epidemiological, and 
economic evidence in health care and qualitative insight in making better decisions. Model 
based analysis allows to quantify different concepts relating to uncertainty in decision mod¬ 
eling. Therefore effective simulation methods like the MOR can contribute in for more pre¬ 
cise decisions and better treatment of diseases. Beside the MOR, the two important fields 
of the modern scientific computing are uncertainity quantification and optimal control. The 
FHN equation is an ideal model in this respect to recover the uncertainities with respect to 
parameters and input variables in the context of neuron modeling. Also optimal treatment of 
the diseases can be designed using reduced order simulation effectively. Optimization and 
uncertainty quantification techniques combined with MOR can improve model predictions, 
to evaluate monitoring schemes and apply better therapy in the medical practice. 
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